Create a Base Map

A male Kentucky Arrow Darter collected from Wild Dog Creek.

First, lets create a kentucky state basemap

world_map <- map_data("world")
county <- map_data("county")
state <- map_data("state")
#apsu_point <- data.frame("x" = -87.353069, "y" = 36.533654)
ky <- county %>% 
  filter(region=="kentucky")

ggplot() + geom_polygon(data = state, aes(x=long, y = lat, group = group),
                        fill = "white", color="black") + 
           geom_polygon(data = ky, aes(x=long, y = lat, group = group),
                        fill = "gray", color="black") + 
  xlab("Longitude") + ylab("Latitude") + ggtitle("Kentucky")

Amazing. Now we have the state of Kentucky, lets add some more data.

kad_co_streams<- sf::read_sf(("/Users/alexisculley/Desktop/aculley_KAD_master_folder/KAD_base_maps/kadstreamdata/kad_stream_data.shp"))
kad_watershed<- sf::read_sf(("/Users/alexisculley/Desktop/aculley_KAD_master_folder/KAD_base_maps/HUC8_CONUS/HUC8_US.shp"))
kad_local <- read.csv("/Users/alexisculley/Desktop/aculley_KAD_master_folder/KAD_base_maps/kad_localities.csv")
sur_mine <-  sf::read_sf(("/Users/alexisculley/Desktop/aculley_KAD_master_folder/KAD_base_maps/Shapefile/5072/2019_activeMining_5072.shp"))

we just added stream, watershed, KAD locality, and surface mining data.

Now lets practice getting data straight from the internet! Lets get a kml file from github that I added.

add .kml file from the web:

local.kml <- readOGR("https://raw.githubusercontent.com/avculley/more_maps/main/kadlocal.kml")
## OGR data source with driver: KML 
## Source: "https://raw.githubusercontent.com/avculley/more_maps/main/kadlocal.kml", layer: "KADlocalities.csv Events"
## with 12 features
## It has 2 fields

Lets make this file usable

local_points <- cbind(local.kml@data, local.kml@coords)
local_points[2] <- NULL
local_points[4] <- NULL
colnames(local_points) <- c("name","x","y")

x = ggplot()+ 
  geom_polygon(data = state, aes(x=long, y = lat, group = group), fill = "white", color="black") + 
  geom_polygon(data = ky, aes(x=long, y = lat, group = group), fill = "gray", color="black") +
    geom_point(data=local_points, aes(x = x, y = y, color = name), size = 4, alpha = 0.8) +
  coord_fixed(xlim = c(-90, -82),  ylim = c(36.2, 39.2), ratio = 1.2)

x

yay they worked!

Now we need to graph these to see what they look like. We added a datset including watershed outlines, localities, and streams.

ggplot(kad_co_streams)+geom_sf()

kad_co_rivers <- subset(kad_co_streams, NAME == "South Fork Kentucky River" | NAME == "Middle Fork Kentucky River" | NAME == "North Fork Kentucky River" | NAME == "Kentucky River")

KAD streams looks great! But.. Hold up, we cant just plot this because it will take 7 years since this dataset contains the entire United States. Lets filter the watershed data to only include the areas of Kentucky we care about. You can do this is a lot of ways, but I am just going to subset the dataset by ID’s I know from the dataset.

HUC8_KAD <- subset(kad_watershed, HUC8 == "05100201" | HUC8 == "05100202" | HUC8 == "05100203" | HUC8 == "05100204")
ggplot(HUC8_KAD)+geom_sf()

Now that we have data on the scale we need, we have proof they are actually into R.

Let’s plot them on top of each other. First, we need to set the crs for the datasets so they match.

st_crs(kad_co_streams)
## Coordinate Reference System: NA
kad_co_streams_crs = st_set_crs(kad_co_streams, 4326)
st_crs(kad_co_rivers)
## Coordinate Reference System: NA
kad_co_rivers_crs = st_set_crs(kad_co_rivers, 4326)
st_crs(sur_mine)
## Coordinate Reference System:
##   User input: NAD83(NSRS2007) / Conus Albers 
##   wkt:
## PROJCRS["NAD83(NSRS2007) / Conus Albers",
##     BASEGEOGCRS["NAD83(NSRS2007)",
##         DATUM["NAD83 (National Spatial Reference System 2007)",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4759]],
##     CONVERSION["Conus Albers",
##         METHOD["Albers Equal Area",
##             ID["EPSG",9822]],
##         PARAMETER["Latitude of false origin",23,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",-96,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",29.5,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",45.5,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Data analysis and small scale data presentation for contiguous lower 48 states."],
##         AREA["United States (USA) - CONUS onshore - Alabama; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming."],
##         BBOX[24.41,-124.79,49.38,-66.91]],
##     ID["EPSG",5072]]
sur_mine_crs = st_set_crs(sur_mine, 4326)
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
## that
ggplot() + 
  geom_sf(data = HUC8_KAD, color = "black") +
  geom_sf(data = kad_co_streams_crs, color = "blue") +
  ggtitle("Kentucky Arrow Darter Historical Range") + 
  coord_sf()

Now that we can plot them on top of each other, lets try and make it look less awful. Lets clip the stream data by the watershed polygons.

kad_watershed_streams = st_intersection(kad_co_streams_crs, st_make_valid(HUC8_KAD))
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
kad_rivers_streams = st_intersection(kad_co_rivers_crs, st_make_valid(HUC8_KAD))
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
test1 = ggplot() + 
  geom_sf(data = HUC8_KAD, color = "black") +
  geom_sf(data = kad_watershed_streams, color = "blue") +
  ggtitle("Kentucky Arrow Darter Historical Range") + 
  coord_sf()

test1

Wow! Now we have all the streams in the watersheds the Kentucky Arrow Darter are found in. This doesn’t tell us a lot though does it?

localmap = ggplot() + 
  geom_sf(data = HUC8_KAD, color = "#CC79A7", fill = "white") +
  geom_sf(data = kad_watershed_streams, color = "#56B4E9") +
  geom_sf(data = kad_rivers_streams, color = "#0072B2") +
      geom_point(data=local_points, aes(x = x, y = y, color = name), color="black", size = 4, alpha = .8) +
      geom_text(data=kad_local,aes(x=long,y=lat,label=name), color="white", vjust=1.6, hjust=-.11, size=2.5, fontface="bold") +
  geom_text(data=kad_local,aes(x=long,y=lat,label=name), color="black", vjust=1.5, hjust=-.1, size=2.5, fontface="bold") +
  labs(x="Longitude",y="Latitude", title="Kentucky Arrow Darter Localities", fill = "No. of Stations") +
  north(kad_watershed_streams, location = "bottomleft", scale = 0.05, symbol = 12) +
  theme(legend.position = "bottom")+
    annotation_scale(location = "bl", width_hint = 0.5) +
  theme(legend.position = "right") +
  coord_sf()
  
#other north arrows
#annotation_north_arrow(location = "br", which_north = "true", 
       # pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
       # style = north_arrow_fancy_orienteering)
#other kad locality file
#geom_point(data = kad_local, aes(x = long, y = lat), color="black", size = 3, alpha = 0.8) +

localmap

Now time to make the inset map.

state <- map_data("state")
county <- map_data("county")
#apsu_point <- data.frame("x" = -87.353069, "y" = 36.533654)
ky <- county %>% 
  filter(region=="kentucky")

ken_ref_map = ggplot() + 
  geom_polygon(data = state, aes(x=long, y = lat, group = group),
                        fill = "white", color="black") + 
           geom_polygon(data = ky, aes(x=long, y = lat, group = group),
                        fill = "gray", color="black") + 
  geom_sf(data = HUC8_KAD, color = "#CC79A7") +
  coord_sf(xlim = c(-90, -82),  ylim = c(36.5, 39)) + 
  theme(panel.background = element_rect(fill = "lightblue"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
        axis.line=element_blank(), axis.text.x=element_blank(), axis.text.y=element_blank(),axis.ticks=element_blank(), 
        axis.title.x=element_blank(), axis.title.y=element_blank())
ken_ref_map

Now lets combine our kad locality map with our inset map!

gg_inset_map1 = ggdraw() +
  draw_plot(localmap) +
  draw_plot(ken_ref_map, x = 0.05, y = 0.65, width = 0.3, height = 0.3)

gg_inset_map1

leaflet(kad_local) %>% 
  addTiles()%>%  
  addCircleMarkers(data = kad_local,
             lat = ~lat, lng = ~long, label = kad_local$name,
                   weight = 2,
                   color = "grey",
                   fillColor = kad_local$name,
                   fillOpacity = 0.7,
             stroke = FALSE, opacity = 0.3)%>%
  addWMSTiles("https://basemap.nationalmap.gov/arcgis/services/USGSTopo/MapServer/WmsServer", layers = "0") %>%
  addMiniMap(zoomLevelOffset = -4) %>%
  addScaleBar()

###it is time to create resistance maps!

Lets begin by plotting surface mining data on the existing KAD map we have made.

ggplot() + 
  geom_sf(data = HUC8_KAD, color = "#CC79A7", fill = "white") +
  geom_sf(data = kad_watershed_streams, color = "#56B4E9") +
  geom_sf(data = sur_mine, color = "red") +
  geom_sf(data = kad_rivers_streams, color = "#0072B2") +
      geom_point(data=local_points, aes(x = x, y = y, color = name), color="black", size = 4, alpha = .8) +
      geom_text(data=kad_local,aes(x=long,y=lat,label=name), color="white", vjust=1.6, hjust=-.11, size=2.5, fontface="bold") +
  geom_text(data=kad_local,aes(x=long,y=lat,label=name), color="black", vjust=1.5, hjust=-.1, size=2.5, fontface="bold") +
  labs(x="Longitude",y="Latitude", title="Kentucky Arrow Darter Localities", fill = "No. of Stations") +
  north(kad_watershed_streams, location = "bottomleft", scale = 0.05, symbol = 12) +
  theme(legend.position = "bottom")+
    annotation_scale(location = "bl", width_hint = 0.5) +
  theme(legend.position = "right") +
  coord_sf()

Yikes! We have way too much mining data. We really only need what is included in the KAD watershed. Lets clip the 2019 surface mining data by the watershed polygons.